DEM 7093 - Some Raster analysis using R

In this analysis, we use R to read in LandSat8 data, create a vegitation index and calculate zonal statistics of a raster. The data come from the 2006 National Land Cover Database.

First we relcassify the NLCD data into two classes, based on the value of the raster. In this case, I want to classify the pixels as to whether they are developed or undeveloped.

We then use a census tract polygon layer to calculate the proportion of each tract’s land area that is developed vs undeveloped.

Raster read in

library(rgdal)

## Loading required package: sp

## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-28, (SVN revision 1158)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Program Files/R/R-4.1.2/library/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Program Files/R/R-4.1.2/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-6
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Overwritten PROJ_LIB was C:/Program Files/R/R-4.1.2/library/rgdal/proj

library(raster)
library(gdalUtils)
nlcd<-raster("C:/Users/ozd504//OneDrive - University of Texas at San Antonio/classes/dem7093/dem7093_22/data/nlcd16_48_lc/nlcd_2016_land_cover_48_20190424.img.vat.img")

nlcdproj<-gdalwarp(srcfile ="C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/nlcd16_48_lc/nlcd_2016_land_cover_change_index_48_20190424.img.vat.img",
                    dstfile = "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/outputnlcd.tif",
                    t_srs = nlcd,
                    output_Raster = T)

## Warning in system(cmd, intern = TRUE): running command '"C:\Program
## Files\QGIS 3.16.16\bin\gdalwarp.exe" -t_srs "raster(ncols=44253,
## nrows=41260, xmn=-1075365, xmx=252225, ymn=309855, ymx=1547655,
## crs='+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0
## +y_0=0 +datum=WGS84 +units=m +no_defs')" -of "GTiff" "C:/Users/ozd504/
## OneDrive - University of Texas at San Antonio/gis_classwork/nlcd16_48_lc/
## nlcd_2016_land_cover_change_index_48_20190424.img.vat.img" "C:/Users/ozd504/
## OneDrive - University of Texas at San Antonio/gis_classwork/outputnlcd.tif"' had
## status 1

raster::plot(nlcdproj,
             main="Original Raster image")

Clip raster to polygon

library(sf)

## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE

## 
## Attaching package: 'sf'

## The following object is masked from 'package:gdalUtils':
## 
##     gdal_rasterize

library(tigris)

## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.

library(terra)

## terra 1.5.21

## 
## Attaching package: 'terra'

## The following object is masked from 'package:rgdal':
## 
##     project

options(tigris_class = "sf")
counties<-tracts(state="TX",
                 county = c("Bexar", "Comal"),
                 cb=T,
                 year=2016)

counties<-st_transform(counties, crs=2278)

# doesn't work

est<-terra::crop(nlcdproj, counties)
plot(est)

est2<-terra::mask(est, counties)
plot(est2)

# est<-mask(x=nlcdproj,mask=as_Spatial(counties))
# est2<-crop(x=nlcdproj, y=extent(as_Spatial(counties)))
# plot(est2)

Raster reclassification:

m<-c(1, 12, 0,
     13, 24, 1,
     25, 95, 0)

rclmat<-matrix(m,
               ncol=3,
               byrow = T)

nlcdreclass<-raster::reclassify(x = est2,
                                rcl = rclmat)

plot(nlcdreclass,
     main="Reclassified Raster image \n 1=developed, 0 = undeveloped",
     col= gray(100:0 / 100))

Load the census tract polygons

library(tigris)
satract<-tracts(state="48",
                county="029",
                year=2018)

library(sf)
satract<-st_as_sf(satract)
satract<-st_transform(satract,
                      crs=  2278)

satract<-as(satract, "Spatial")
summary(satract)

Calculate zonal statistics

Now we find the mean of the binary reclassified data, to estimate the percent of each tract’s area that is developed.

First we reproject the raster to the same coordinate system as the census data. I get the projection string from ^^^ above.

nlcdreclass<-projectRaster(nlcdreclass,
                           crs="+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0")

satract$pct_developed<-as.numeric(extract(nlcdreclass, satract, fun=mean)) #this takes a minute or two

## Warning in .local(x, y, ...): Transforming SpatialPolygons to the crs of the
## Raster

satract<-st_as_sf(satract)

library(ggplot2)

## 
## Attaching package: 'ggplot2'

## The following object is masked from 'package:terra':
## 
##     arrow

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:terra':
## 
##     intersect, src, union

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

satract<-satract%>%
  filter(complete.cases(pct_developed))

tractlines<-st_boundary(satract)

library(tmap)
tm_shape(satract)+
  tm_polygons("pct_developed",style="quantile", n=5,legend.hist=T )+
  tm_shape(tractlines)+
  tm_lines(col="black")+
  tm_format("World",
            title="San Antonio Developed Land Area",
            legend.outside=TRUE )

I used this method in a paper I wrote in 2013, although in that paper I did the diversity of the land cover. That could be done like:

library(vegan)

## Warning: package 'vegan' was built under R version 4.1.3

## Loading required package: permute

## Warning: package 'permute' was built under R version 4.1.3

## Loading required package: lattice

## This is vegan 2.5-7

shannonVegan <- function(x, ...) {
    diversity(table(x), index="shannon")
}

#shanVegOut <- focal(nlcdproj, fun=shannon, pad=T)  

satract$var_land<-as.numeric(extract(nlcdproj,
                                     satract,
                                     fun=shannonVegan)) #this takes a minute or two

shannon<-tm_shape(satract)+
  tm_polygons("var_land",
              style="quantile",
              n=5, 
              legend.hist=T )+
  tm_shape(tractlines)+
  tm_lines(col="black")+
  tm_format("World",
            title="San Antonio Shannon Diversity of Land Cover",
            legend.outside=TRUE )

shannon

Exporting images for posters/presentations

I recommend using the ggsave() function to export images made by ggplot for use in other documents. For tmap, you can use tmap_save()

#ggsave(m1, path = "", filename = "sa_pct_devel.png", units = "in", width=8, height = 10)
tmap_save(shannon,
          filename = "~/OneDrive - University of Texas at San Antonio/classes/gis_classwork/shannon_tract.png",
          width = 8,
          height = 10,
          units = "in")

And you can insert that image into a doc or presentation.

Landsat imagery

Good simple description of landsat here

library(raster)
library(rgdal)
library(rgeos)

## rgeos version: 0.5-9, (SVN revision 684)
##  GEOS runtime version: 3.9.1-CAPI-1.14.2 
##  Please note that rgeos will be retired by the end of 2023,
## plan transition to sf functions using GEOS at your earliest convenience.
##  GEOS using OverlayNG
##  Linking to sp version: 1.4-6 
##  Polygon checking: TRUE

library(RColorBrewer)


# turn off factors
options(stringsAsFactors = FALSE)
landsatstack<- list.files("C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/", 
  pattern=glob2rx("*B*.TIF$"), 
  full.names = T)

landsatstack

##  [1] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B1.TIF" 
##  [2] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B10.TIF"
##  [3] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B11.TIF"
##  [4] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B2.TIF" 
##  [5] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B3.TIF" 
##  [6] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B4.TIF" 
##  [7] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B5.TIF" 
##  [8] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B6.TIF" 
##  [9] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B7.TIF" 
## [10] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B8.TIF" 
## [11] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_B9.TIF" 
## [12] "C:/Users/ozd504/OneDrive - University of Texas at San Antonio/gis_classwork/LC08_L1TP_027040_20200311_20200325_01_T1/LC08_L1TP_027040_20200311_20200325_01_T1_BQA.TIF"

lsb2<-raster(landsatstack[2])

plot(lsb2,col = gray(0:100 / 100))  

  rm(lsb2)      ; gc()

##            used  (Mb) gc trigger  (Mb)  max used   (Mb)
## Ncells  4349763 232.4    8465130 452.1   8465130  452.1
## Vcells 26157396 199.6  101976122 778.1 451409981 3444.0

Make raster stack of bands 1 - 7

lsstack<-stack(landsatstack[c(1,4:9)])

lsbrick<-brick(lsstack)
names(lsbrick)

## [1] "LC08_L1TP_027040_20200311_20200325_01_T1_B1"
## [2] "LC08_L1TP_027040_20200311_20200325_01_T1_B2"
## [3] "LC08_L1TP_027040_20200311_20200325_01_T1_B3"
## [4] "LC08_L1TP_027040_20200311_20200325_01_T1_B4"
## [5] "LC08_L1TP_027040_20200311_20200325_01_T1_B5"
## [6] "LC08_L1TP_027040_20200311_20200325_01_T1_B6"
## [7] "LC08_L1TP_027040_20200311_20200325_01_T1_B7"

names(lsbrick) <- gsub(pattern = "LC08_L1TP_027040_20200311_20200325_01_T1_",
                       replacement = "",
                       names(lsbrick))
names(lsbrick)

## [1] "B1" "B2" "B3" "B4" "B5" "B6" "B7"

Normalized Density Vegitation Index - NDVI

ndvi<-(lsbrick[[5]] - lsbrick[[4]]) / (lsbrick[[5]] + lsbrick[[4]])

# plot(ndvi,
#      main = "NDVI around San Antonio",
#      axes = FALSE, box = FALSE)       
library(tmaptools)
tm_shape(ndvi)+
  tm_raster(style="quantile",
            n=5,
            legend.hist=T,
            palette =get_brewer_pal("RdYlGn", n = 5) )+
  tm_format("World",
            title="NDVI - San  Antonio Region",
            legend.outside=TRUE )+
tm_shape(satract)+
  tm_polygons(alpha = .7)

## stars object downsampled to 991 by 1010 cells. See tm_shape manual (argument raster.downsample)

## Variable(s) "NA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

hist(ndvi)

library(mapview)
mapview(nlcdreclass)

## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 9494355 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  9494355 '

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition

writeRaster(x= ndvi, 
    filename = "C:/Users/ozd504//OneDrive - University of Texas at San Antonio/gis_classwork/ndvi1_22.tif", format="GTiff", datatype="INT2S", overwrite=T)

Landsat plots

plotRGB(lsbrick,
r = 6, g = 5, b = 2,
stretch = "lin",
axes = TRUE,
main = "RGB composite image\n showing Agricultural land - Landsat Bands 6, 5, 2")
box(col = "white")

plotRGB(lsbrick,
r = 5, g = 4, b = 3,
stretch = "lin",
axes = TRUE,
main = "Color infrared composite image\n Landsat Bands 5, 4, 3")
box(col = "white")